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We investigate states of enhanced activity in a biological neuronal network composed of pulse-coupled oscillators. 
The synaptic couplings between the neurons are dynamic, modeling spike time dependent plasticity. The network 
exhibits statistical characteristics which recently have been identified in an analysis of epileptic seizures [Osorio 
et al., Phys. Rev. E82, 021919(2010)] based on analogies to the onset of earth quakes. 
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I. INTRODUCTION 

The origin of epileptic seizures is still a mystery (see e.g. HI). 
The epileptic brain exhibits spatio-temporal complexity, both 
during seizures as well as in the inter-seizure intervals. In 
recent years, the dynamical properties of inter-seizure states 
have been in the focus of research with the aim to formulate 
suitable measures for the prediction of seizure onset (2). 

A major goal of theoretical seizure modeling is the develop- 
ment of models of neuronal activity, which are able to repro- 
duce the observations and allow for a theoretical interpretation 
of the spatio-temporal behaviour of seizure states. A successful 
model will shed considerable light on the questions whether 
seizures are predictable and what measures can be rewarding. 
The basic issues of theoretical seizure modeling are to inves- 
tigate how epileptic seizures emerge, spread, and terminate 
in such models. Dynamical systems theory points to the im- 
portance of synchronization, however, it seems to be rather 
unclear whether the synchronization of neuronal activity is 
a by-product, whether it induces or whether it terminates a 
seizure 0H]]. 

Recently, Osorio et al. opened up a new perspective for 
our understanding of epilepsy assessing statistical similarities 
between epileptic seizures and the onset of earthquakes. Based 
on a comparative study of seizure and earth quake data, the 
authors showed ample evidence that the Gutenberg-Richter law, 
the Omori law and the statistics of interevent times formulated 
for earth quake data can be detected in epileptic seizure data as 
well. 

The purpose of the present article is to develop a neural net- 
work model, which is able to reproduce the main characteristics 
of the analysis of Osorio et al. [5]. The developed model is 
based on a generalization of Haken's Lighthouse model (for a 
review we refer the reader to [6]), which is augmented by a suit- 
able learning algorithm taking into account the phenomenon of 
spike-timing-dependend plasticity (STDP), following closely 
the work of Chen et al. (73 ■ The model presented in this article 
can be viewed as a model for epileptic seizures similar to the 
Burridge-Knoppoff model [8 1 and its variants for earth quakes. 
Our extension of the Lighthouse model allows for the inclu- 
sion of a spontaneous restructuring of the network architecture, 
making use of the idea that dynamical seizures of the brain can 
be learned or dislearned [9|, 1 10 1 . 

It is well-known that adaptive networks may lead to the emer- 
gence of critical states in neural network models based on the 
phenomenon of self-organized criticality (SOC), as described 
by P. Bak ifTTl . Experimentally, Beggs and Plenz lfl2l showed 
convincing evidence that the dynamics of neuronal popula- 
tions could actually exhibit critical dynamics by an explicit 
investigation of events called neuronal avalanches. Neuronal 
avalanches can be found in regions of control parameter space 
located in between a regime of coherent, wavelike neuronal 
activity and a regime, where neurons behave asynchronous 



and spike incoherently. The existence of SOC in an integrate 
and fire model including dynamical synapses has been shown 
by A. Levina et al. lfT3l [T4l . The dynamics of their model 
exhibits the occurence of critical avalanches and the distribu- 
tion of the size of the occurred avalanches showed a power 
law behaviour. In the case of static synapses, the coupling 
strengths have to be fine-tuned in order to generate a critical 
state, whereas dynamical synapses renders the system critical. 
Similar results were obtained by L. de Arcangelis et al. 1731 . 
C. Meisel et al. [16 | presented a study of neural networks 
including spike-time-dependent synaptic plasticity, extending 
the earlier work of Bornholdt and Rohlf IfTTl . They showed 
that due to the inclusion of this learning mechanism the net- 
work self-organizes into a critical state. Millman et al. 1T81 
demonstrated the existence of SOC in nonconservative models 
of networks of leaky integrate and fire neurons with short-term 
synaptic depression. They demonstrated the existence of two 
states, up and down states and demonstrated that up states are 
critical, whereas down states are subcritical. 

With respect to epileptic seizures it is not clear whether SOC 
plays a dominant role. Recently, it has been discussed that 
adaptive self-organized criticality fails during epileptic seizure 
attacks ITU 

The present article is outlined as follows. In the first section, 
we review Haken's Lighthouse model and discuss its exten- 
sion taking into account a learning mechanism based on the 
mechanism of spike-timing-dependent plasticity. In the second 
section, we investigate in detail states of enhanced activity. We 
reproduce features similar to the ones presented in lfl4l l20l . 
Then, we focus on similarities with the investigations of exper- 
imental recordings of seizure states discussed in 0. A central 
issue is the question why states of enhanced activity emerge 
and terminate in the used model. 



II. MODEL AND IMPLEMENTATION 

A. Lighthouse Model 

The Lighthouse model has been introduced by Herman 
Haken and is described in great detail in his monograph |5). It 
models brain tissue as a neuronal network composed of pulse 
coupled oscillators. Each neuronal oscillator, labeled by the in- 
dex m, is characterized by its phase <f) m , a 2n- periodic variable, 
and a dendritic current tf/ m . 

The phase describes the action potential of a neuron. If the 
phase reaches the threshold 2n, a neuron will spike. The action 
potential S k generated by neuron k is obtained as a sum over 
delta functions: 
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The temporal evolution of the phase is determined by the equa- 
tion 



(f> m (i) = S 



P(ext,m) (0,® 

, k 
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where £<(Ek c mk t /'k(t) + P( e xt,m)(t),®) is a sigmoidal function of 
the sum of all dendritic currents t/^ of neurons k mutiplied by 
a weighting factor c m k and an external signal P( e xt,m)(t)- The 
factors c m £ form the increment matrix c. In the following, we 
use a diagonal increment matrix c m k = c m 5 m k. The sigmoidal 
function E(X,@) is taken as the Naka-Rushton relation and 
describes the response of a neuron to an input current. 
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For high input values X the output saturates at a maximal firing 
rate v. If the input signal is below a certain threshold value 
®, the output will tend to zero. As outlined in [6 1 this models 
the all-or-none-behaviour of the axon hill, which generates the 
action potential. Due to the saturation of the Naka-Rushton 
relation for large input values, a refractory period of the neurons 
is included in the system. The constant M is related to the slope 
of the function. In our treatment we have used M — 3 and a 
maximal firing rate v of 1 . 

When a neuron generates an action potential, this will be 
transfered to all connected neurons. The dendritic current if/ m (t) 
of neuron m will change due to the input from neuron k. It 
obeys the following differential equation 



>frm(t) = ^ a mk S k {(p k (t)) - yifr m (f). 

k 



(4) 



The synaptic weights are denoted by the matrix elements a mk 
and describe the strength of the connection from neuron k to 
neuron m. The synaptic weights characterize the connections of 
the neural network. In the following self-coupling is neglected, 



inter-spike intervals tj+\ — tj. The coupling weight aij decreases 
in the first (acausal) case. A similar reasoning can be used for 
the reversed cased. It is convenient to introduce two functions 
<Tj(t), denoted as A;(f) and Bdf), for each neuron j. 

The variable A m (t) associated to neuron m is related to 
strengthening of the connection a^m if the firing with respect 
to neuron k is causal, i.e. if the product A m (t)Sk(t) is differ- 
ent from zero. The variable Bkit) induces weakening of the 
coupling strength a km if the firing with respect to neuron m is 
acausal, i.e if the product Bk(t)S,„(t) is different from zero. The 
dynamics of the coupling weight is then given by the evolution 
equation 



hk m - AA m S k - ra km B k S 



(5) 



The first term of equation Q is related to strengthening due to 
causal firing and the amount is proportional to the potentiation 
constant A. The second term causes depression of the coupling 
weight. It is proportional to the depression constant r and to 
the coupling weight ak,„ itself. 

In order to determine the functions A„,(f) and B m (t) dynami- 
cally we closely follow the work of Chen and Jasnow [7 1. They 
used the following evolution equations, where cr m (t) now de- 
notes A m and B m . The quantities A m and B m can be viewed as 
concentrations of chemical species, which are generated by a 
spiking event of the neuron m. 



& m =u (T (l-cr m -I )S m ; cr-A,B (6) 

To- 

Furthermore, the dynamics of is given by 

cr I a 
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i.e. a„ 



:0. 



The spiking of a neuron, encoded in S m , leads to an increase 
of the concentrations A m , B m . Furthermore, the concentrations 
decay with time constants tv. For details we refer the reader to 
0. 



B. Spike-Timing-Dependent Plasticity 

In order to include learning, we combine the Lighthouse 
model with the mechanism of spike time dependent plasticity. 
For a historical review we refer the reader to ETl . This mecha- 
nism is based on the Hebbian learning rule. We will rely on the 
formulations of C. Chen and D. Jasnow Q and van Rossum et 
al. El . The underlying mechanism of STDP is as follows: 
The coupling weight ak m will be strengthened due to causal 
firing and weakened in the case of acausal firing. Causal firing 
means that neuron m spikes in advance of neuron k. Acausal 
firing denotes the reverse case. The reverse coupling weight 
a,„k will be weakened if a km is strengthened and vice versa. 
Hence, the connections between two neurons become directed 
and dynamical. 

In order to formulate a mathematical model of spike time 
dependent plasticity one has to design a mechanism, which is 
able to distinguish causal from acausal spiking. Spiking events 
of neuron i are encoded in the spike train S ; = 6{t - t"). To 
the n-th firing event of neuron j at time fl we connect a function 
u-,".) 

<Tj(t) = <xo0(f - tj)e~ 7 <r . The product cr ; (f) 5,(f) is zero if 

t" < t", whereas it is nonzero in the reverse case. We have to 

' j 

assume that the decay time T a is smaller than the characteristic 



in. STATES OF ENHANCED ACTIVITY 

The aim of this article is to find states in the extended Light- 
house model, which exhibit features of epileptic seizures. There 
seems to be no commonly accepted definition of a seizure in 
literature. However, as a main feature one may consider the 
fact that all neurons in a region with epileptic spiking are active 
and spike with a high frequency. Synchronization in the sense 
of dynamical systems theory seems to be an acompanying phe- 
nomenon, as pointed out by R. S. Fisher et al l23l . However, 
whether synchronization is the origin of a seizure or plays a 
major role in the terminiation, as emphasized by |H][24], is yet 
unclear. 

Therefore, in our analysis of the generalized Lighthouse 
model, we use the notion of states of enhanced activity, which 
could be interpreted as an epileptic seizure. A state of enhanced 
activity is defined by a strong mean dendritic current 

V(t) = ^^Mt) (8) 

i 

A plot of the time signal of this mean current is exhibited in fig 
Q. It is evident that periods of strong mean current and weak 
mean current interchange. In order to explicitly determine 
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Figure 1 : (a) Phase (p of a single neuron is shown as a function of 
time, (b) The dendritic current if/ of the same neuron is plotted in time. 
Figures (a) and (b) demonstrate the existence of time-intervals, during 
which the neuron is at rest and intervals when it is higly active. 



a state of enhanced activity we introduce a threshold value 
for the mean current. A closer inspection of the time signals 
demonstrates that during the state of enhanced activity the 
spiking frequency of most of the neurons are in the saturated 
regime, i.e. almost all neurons fire with the highest attainable 
frequency v. This is demonstrated in fig Q, in which the mean 
spiking frequency < ip > 

<<*(0>=^J>(0 (9) 

i 

of the system is shown. 



A. Initiation and Termination of States of Enhanced Activity 

States of enhanced activity emerge in the neuronal generalized 
Lighthouse model in the following situations: The coupling 
weights are initialized randomly which implies that all to all 
coupling exists. One special neuron is singled out receiving an 
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Figure 2: (a) Mean dendritic current *P(f) as a function of time, (b) 
Mean spiking frequency < tp > as a function of time. 

external signal P( e xt,ni) consisting of a sequence of delta peaks: 
P(ex,,m) = J]A Pe J(t - nAT) (10) 

n 

The amplitude A Pext of the peaks is chosen in such a way that 
the externally excited neuron always gains the maximal spiking 
frequency v. The period AT is choosen to be 2n/v. 

The external signal mimics an external stimulation of the 
network. It is well-known that epileptic seizures can be initiated 
by periodic external stimuli like light flashes. 
As one can see in figure ([TJ, time-intervals occur where all 
neurons emit nonvanishing dendritic currents connected with a 
rapidly changing phase implying a strong spiking behaviour of 
the neurons. These time intervals of strong spiking activity are 
interupted by intervals where no neuron, except the externally 
driven one, is active. 

States of enhanced activity are initiated by the externally 
driven neuron. Due to the external forcing this neuron spikes 
periodically. As a consequence, all neurons receive pulses, 
whose strength is proportional to the coupling weight. The 
phases of the neurons pile up until the phase of one of the neu- 
rons reaches 2n. Then, this neuron firing leads to an increase 
of the phase of the other connected neurons. If their phases are 
close to 2n, they eventually will spike and a chain reaction is 
initiated, which results in a state of enhanced activity. 
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If the amplitude A Pea is lowered, the seizure onset time 
increases. Below a certain threshold, states of enhanced activity 
are not observed. 

Termination of a states of enhanced activity is clearly due 
to a change of the coupling weights a^. During states of 
enhanced activity the causal nature of the firing of neurons, 
which is responsible for the built up of a network structure 
with strong couplings between causally firing neurons, change 
to acausal firing. This eventually leads to a weakening of the 
coupling weights and, hence, to the destruction of the state of 
enhanced activity. 

Let us exemplify this mechanism considering two neurons i 
and j in more details. Assume that neuron j fires before neuron 
i. As a result the coupling weight will increase, whereas a» 
decreases. Hence, we investigate the following dynamics: 



4>i = 3(^,0) 

h = -yi!', + (i,jS j((bj) 

fa = 3(^,0) 



Numerical parameters: 


Naka-Rushton relation: 


Time step dt = 0.01 


maximal spiking rate v: 1.0 




i nresnoiu a. iu.u 




oLCCpncso IV1 . J 


External signal: 


Neuronal network: 


Amplitude A Pa : 10.0 


Number Neurons N: 50 


Forcing Periodic T Pea : 1.0 


Damping constant y: 0.7 




Enhancement constant c m : 5.0 


Learning parameters: 


Potentiation constant A: 1 .0 Depression constant r. 1 .0 


Potentiation width r^: 0.2 


Depression width t#: 0.2 


Recovery rate r rcr : 10.0 


Fatigue rate T/o-: 10.0 


Chemical fraction u^: 0.9 


Chemical fraction ug: 0.9 



(11) 



Table I: Numerical coefficients used in the simulations, are shown. 



We have neglected the coupling ajj w 0, taking into account 
that the neuron j spikes before neuron i. 

The solution for the dendritic current \j/i(i) is given in terms 
of the spike-times T" of neuron j, which are determined by the 
spiking condition © ; (T") = 2n: 



(MX) = 2* 



-70- 



-T") 

1 a; 



(T») 



(12) 



If we consider a step function E(ip, 0), which is zero if < 

and y, if tfr > 0, the phase fy(t) changes only, if > 0. This is 

the case if a,- ,<T?) >0. 
j j 

Due to the acausal situation of neuron j, i.e. due to ap « 
the dendrictic current of neuron j decays exponentially, and the 
spiking frequency decreases 



<f>i = <t>. 



; (0) + dt'E 



0/(0) + vln 



1 | ( p®l }M 



-Myt^ 



(13) 



Eventually, this leads to a stationary phase <pj(Q) + 

vln^(^^-) M ). However, for small values of the damping con- 
stant y, the dendritic current of neuron i remains at a finite 
value, and, in turn its spiking frequency remains finite, pro- 
vided the condition Y„ e~ yi1 ~ T '] ) aij (T n ) > is fullfilled. As 
a consequence, the phase (£>,- increases more rapidly than the 
phase of neuron j: The phase of neuron i overtakes the phase of 
neuron i rendering the connection ay acausal. This transition 
from causal to acausal firing leads to a decrease of the coupling 
weight flji and an increase of the coupling a». 

This observation can be considered to be consistent with the 
experimental findings of Schindler et al. [24 1 and Lehnertz et 
al. |H|. These authors state that for a certain kind of epileptic 
seizures synchronization and instability occur simultaneously. 
They argue that the synchronization towards the end of the 
seizure may be even a self -regulatory mechanism to catalyse 
the seizure termination. 

The described behaviour exists for a certain range of values of 
the damping constant y, y l c < y < y". This has been investigated 



by a numerical investigation of the generalized Lighthouse 
model. This range depends on the system size. 

The described mechanism only works, if the damping con- 
stant y of the dendritic current is large enough. Therefore a 
lower critical damping constant y[ exists. On the opposite, the 
chain reaction can not start, if the damping constant exceeds a 
certain critical upper value y". The neurons will start to spike, 
but the dendritic current of the neuron is decaying too fast and 
therefore the phase velocity is not high enough to preserve the 
firing and initialize a state of enhanced activity. 



B. Analysis of Statistical Behaviour 

In the following, we present results of the investigation of the 
statistical behaviour of the states of enhanced activity. In or- 
der to create sufficient statistics long time simulations of the 
generalized Lighthouse model and measurements of the mean 
dendritic currents have been performed. The used numerical 
parameters are shown in table fl). 27 runs have been realized, 
each with 47 • 10 7 timesteps. In total 10,000 events, which 
could be classified as state of enhanced activity , have been 
measured. The duration of a state of enhanced activity, the 
interevent waiting times between two subsequent states of en- 
hanced activity, and the total energy of a state of enhanced 
energy have been calculated, The energy of a state of enhanced 
activity is defined to be proportional to the sum of the squares 
of all dendritic currents: 



E = R 



(14) 



Based on this data, we have calculated the probability den- 
sity functions (PDF) of the energy, the event duration and the 
interevent waiting time. Furthermore the Omori- and inverse- 
Omori-Law was determined in the same way as described 
by Osorio et al. [5 |. Furthermore, the conditional expected 
waiting-time was calculated. 

In order to assess the scaling behaviour of the various quanti- 
ties like energy of an state of enhanced activity, duration of an 
state of enhanced activity, or the interevent times between two 
subsequent states of enhanced activity, the scaling coefficients 
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Figure 3: The logarithm of the probability density function /(/) of the 
interevent-intervals of the states of enhanced activity as a function of 
Inl. The inset shows the same data points in /(/)-/ coordinates . A 
power law is fitted to the data. 



of the various quantities are calculated by the method intro- 
duced by Newman 11251 . It is based on a maximum likelihood 
estimation of the scaling coefficient, The estimation is only 
applied in intervals [x m i n ,x max ] of the data in which the power 
law behaviour holds. 



C. Interevent Waiting Time 

As stated by Osorio et al and Saichev et al [5, 26], the proba- 
bility density function of the interevent waiting time for earth- 
quakes is described by a power law with a scaling coefficient 
of —1.1. Osorio et al [5] calculated the PDF for data sets of 
epileptic seizures. They obtained an approximate value of — | 
for the scaling coefficient. 

From our numerical computations of the generalized Light- 
house model we extract a scaling coefficient of - 1 .495 + 0. 149 
for [x m i n - 100,x max = 18000], which perfectly agrees with 
the analysis of Osorio [5 1. The data is exhibited in figure {3). 
It is obvious that the power law fits quite well the measured 
behaviour for small and medium interevent waiting times. For 
larger values of the interevent times the data points are scattered 
around the fitted power law. This could be explained, following 
Osorio et al [5], by the assumption that the probability density 
function exhibits different power law regimes, a possibility 
which has first been pointed out by lf26l . However, it could also 
be attributed to the restricted sample size. We remind the reader 
that the behaviour in the regime of large interevent times is of 
considerable interest with respect to the analysis of extreme 
events in complex systems. 



D. Gutenberg-Richter-Law 

The Gutenberg-Richter-law states that the probability density 
function for the earthquake seismic moment S exhibits a power 
law distribution. Similarly, the analysis for epileptic seizures 
by yields a value close to - 5 /3. Furthermore they showed 
that a power law with the same scaling exponent is obtained 
for the energy distribution of epileptic seizures. The results of 
our calculations, which are summarized in figure (H), support 



Figure 4: The probability density function f(E) of the total energy of 
state of enhanced activity is shown. The coefficient is - 1. 660 ± 0.220 
in the interval [3.7 • 10 15 ,4 • 10 16 ]. 
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Figure 5: The logarithm of the probability density function f(D) of 
the duration as a function of log D. The power law fit yields a value 
for the scaling coefficient of - 1 .659 ± 0. 1 55 in the interval [200, 3783] . 
The inset exhibits f(D) as a function of D. 



this assumption. 

As for the statistics of the interevent times the behaviour of the 
PDF deviates from the powerlaw for high values of E, which 
can be explained by the same reasonings as discussed above. If 
one considers the PDF of the duration of the states of enhanced 
activity, which is plotted in figure Q, one can see that both 
PDFs have a similar scaling exponent. 

Therefore, a direct relationship between energy E and duration 
of the state of enhanced activity should exist. This relationship 
can be seen on the basis of the definition 



£(5)= ^Z^ m(f " )2 = <I>D ( S )- 

t„,t„ 



(15) 



of the energy contained in the state of enhanced activity be- 
tween the time instants t a , t e . If we replace approximately the 
quantity iff m (t n ) 2 by a mean intensity < / > we obtain a direct 
relationship between E and duration D, and, hence, the same 
scaling behaviour. 

Assuming the validity of this relation, one can conclude 
that all states, independent of their duration, have the same 
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Figure 6: The foreshock and aftershock rates for states of enhanced 
activity are shown in a double logarithmic scale. The time axis for 
the foreshock rate is inverted. The Omori-Law is fitted into the data 
points and the fitting-parameters are shown. 



finite size effect and depends on the number of neurons of the 
model. A similar behavior is observed for epileptic seizures as 
well as earthquakes J5]. 
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Figure 7: Averaged waiting time to the next seizure-like-state as a 
function of the time distance to the last event. 



static distribution of the intensity < / > and therefore the same 
distribution of the dendritic current V F. 



IV. SUMMARY AND OUTLOOK 



E. Omori Law and Inverse Omori Law 

Earthquakes are rare events, which are accompanied by the 
occurence of aftershock and foreshock events. The rate n(At) 
of these events follows the so-called Omori- and inverse Omori- 
Law as a function of the time distance At from the main event: 



log(n(Af)) = log(fc) - p log(c + Af) 



(16) 



Both laws are empirical. These laws originally were introduced 
by Omori et al [ 27 1 and expanded by Utsu et al 11281 . 

Osorio et al J5i| pointed out that the Omori- and inverse 
Omori-Law can also be formulated for epileptic seizures. The 
analysis of data generated by the generalized Lighthouse model 
supports the validity of these laws. Our results are presented 
in figure The Omori-Law is fitted to the data. Both, the 
foreshock and aftershock rates can be represented by the same 
fit. This implies that the shock rates are not only described by 
the Omori-Law, but are symmetric. The results are consistent 
with the results for data of epileptic seizures |5 1 and are similar 
to the earthquake data Il5l I29H3T1 . 



F. Expected Waiting Time 

An interesting question is how the waiting time, which is 
the time to the next event, depends on the time that has al- 
ready passed since the last event. The results of our numerical 
simulations are shown in figure ([7]). 

As one can see, the expected waiting time to the next event 
increases as a function of the distance to the last event to a 
maximum. Then the expected waiting time rapidly drops down 
to zero. The intersection with the x-axis marks the maximal 
time that can pass before the next state occurs. This is clearly a 



We combined the Lighthouse model introduced by Haken 
with a STDP-learning algorithm to a model in which states of 
enhanced activity emerge. It has been shown that these states 
share characteristic features reported for time signals of epilep- 
tic seizures. Defining states of enhanced activity via the energy 
E we demonstrated the validity of the Gutenberg-Richter law, 
the Omori and inverse Omori laws and showed similarities in 
the statistics of inter-seizure intervals. We want to emphasize 
that our model yields scaling exponents rather close to the ones 
reported in the work of Osorio et al. [5 1. Therefore, a detailed 
investigation of the generalized Lighthouse model should al- 
low for an analytical assessment of these scaling exponents. 
Furthermore, we have described the mechanism which leads to 
a termination of the states of enhanced activity as a transition 
from causal to acausal firing, connected with an immediate 
restructuring of the network topology due to STDP. Due to 
the closeness to experimental findings it would be interesting 
to explore the phase diagram of this generalized Lighthouse 
model in more details and analyze experimental data, whether 
the mechanism of termination of state of enhanced activity can 
also be detected for epileptic seizures. For the future, we plan 
to perform an extended analysis in order to assess the transi- 
tions and bifurcation scenarios leading to the regime, where 
state of enhanced activity arise, from ordered or disordered 
states adding dynamical details to the schematic phase diagram 
given in Osorio [5] locating seizure states in a control parame- 
ter space spanned by the variables termed interaction strength 
and heterogeneity. Additionally, we plan to extend our analysis 
to larger data samples extending our numerical calculations to 
longer time periods. This will allow us to resolve the PDF's 
of extreme events, along the lines suggested by Sornette f32l . 
In this way we want to contribute to the recent discussions on 
the existence of dragon kings Il32l [33ll and the possibility of 
predictions of extreme events, recently reviewed in 
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